function [ m_ccc,m_css,m_cww,m_cee,m_cnn,m_cne,m_cnw,m_cse, m_csw ] = rrb_hyper( Nx , Ny , cc , cs, cw)
%function [ m_ccc,m_css,m_cww,m_cee,m_cnn] = rrb_index( Nx , Ny , cc , cs, cw)

%Construction of RRB Preconditioner
%   Detailed explanation goes here

m_ccc =cc;
m_css =cs;
m_cww =cw;

m_cnn = zeros(Nx+3,Ny+3);
m_cee = zeros(Nx+3,Ny+3);
for jj= 1:Ny+2
    for ii=1:Nx+2
        m_cnn(jj,ii) = - cs (jj+1,ii);
        m_cee (jj,ii) = -cw (jj,ii+1);
    end
end
%generatespy_index(Nx,Ny,m_ccc,m_css,m_cww,m_cee,m_cnn);

m_cne = zeros(Nx+3,Ny+3);
m_cnw = zeros(Nx+3,Ny+3);
m_cse = zeros(Nx+3,Ny+3);
m_csw = zeros(Nx+3,Ny+3);

m_Nx1 = Nx;
m_Nx2 = Ny;

st = 1;
tw = 2;
niv = 0;
while tw < m_Nx2+2
    if (niv > 0)
        fac = 1.0 + (st) * 10.0 / ((m_Nx1 + 2) * (m_Nx2 + 2));
        %fac=1.0;
        for jj = st+1:tw:m_Nx2+3
            for ii = 1:tw:m_Nx1+3
                m_ccc(jj,ii) = m_ccc(jj,ii) * fac +m_cne(jj,ii) + m_cnw(jj,ii) +m_cse(jj,ii) + m_csw(jj,ii);
                m_cne(jj,ii)=0;
                m_cnw(jj,ii)=0;
                m_cse(jj,ii)=0;
                m_csw(jj,ii)=0;
            end
        end
        
        for jj =1:tw:m_Nx2+3
            for ii = st+1:tw:m_Nx1+3
                m_ccc(jj,ii) = m_ccc(jj,ii) * fac +m_cne(jj,ii) + m_cnw(jj,ii) +m_cse(jj,ii) + m_csw(jj,ii);
                m_cne(jj,ii)=0;
                m_cnw(jj,ii)=0;
                m_cse(jj,ii)=0;
                m_csw(jj,ii)=0;
            end
        end
    end
    
    niv=niv+1;
    % Gaussian Elimination: This is done for all R1 and R2 Nodes.
    
    for jj = st+1:tw:m_Nx2+3
        for ii = st+1:tw:m_Nx1+3
            iw = ii - st;
            ie = ii + st;
            js = jj - st;
            jn = jj + st;
            if (jn <= m_Nx2 + 3)
                coeff = m_cnn(jj,ii)/ m_ccc(jn,ii);
                m_cne(jj,ii) = m_cne(jj,ii) - m_cee(jn,ii)*coeff;
                m_cnw(jj,ii) = m_cnw(jj,ii) - m_cww(jn,ii)* coeff;
                m_ccc(jj,ii) = m_ccc(jj,ii) - m_css(jn,ii) *coeff;
                m_cnn(jj,ii) = -m_cnn(jn,ii) *coeff;
            end
            if (js > 1)
                
                m_cse(jj,ii) = m_cse(jj,ii) - m_css(jj,ii) * m_cee(js,ii) / m_ccc(js,ii);
                m_csw(jj,ii) = m_csw(jj,ii) - m_css(jj,ii) * m_cww(js,ii) / m_ccc(js,ii);
                m_ccc(jj,ii) = m_ccc(jj,ii) - m_css(jj,ii) * m_cnn(js,ii) / m_ccc(js,ii);
                m_css(jj,ii) = -m_css(jj,ii) * m_css(js,ii) / m_ccc(js,ii);
            end
            if (ie <= m_Nx1 + 3)
                
                m_cne(jj,ii) = m_cne(jj,ii) - m_cee(jj,ii) * m_cnn(jj,ie) / m_ccc(jj,ie);
                m_cse(jj,ii) = m_cse(jj,ii) - m_cee(jj,ii) * m_css(jj,ie) / m_ccc(jj,ie);
                m_ccc(jj,ii) = m_ccc(jj,ii) - m_cee(jj,ii) * m_cww(jj,ie) / m_ccc(jj,ie);
                m_cee(jj,ii) = -m_cee(jj,ii) * m_cee(jj,ie) / m_ccc(jj,ie);
            end
            if (iw > 1)
                
                m_cnw(jj,ii) = m_cnw(jj,ii) - m_cww(jj,ii) * m_cnn(jj,iw) / m_ccc(jj,iw);
                m_csw(jj,ii) = m_csw(jj,ii) - m_cww(jj,ii) * m_css(jj,iw) / m_ccc(jj,iw);
                m_ccc(jj,ii) = m_ccc(jj,ii) - m_cww(jj,ii) * m_cee(jj,iw) / m_ccc(jj,iw);
                m_cww(jj,ii) = -m_cww(jj,ii) * m_cww(jj,iw) / m_ccc(jj,iw);
            end
        end
    end
    
    for jj = 1:tw:m_Nx2+3
        for ii = 1:tw:m_Nx1+3
            iw = ii - st;
            ie = ii + st;
            js = jj - st;
            jn = jj + st;
            if (jn <= m_Nx2 + 3)
                m_cne(jj,ii) = m_cne(jj,ii) - m_cnn(jj,ii) * m_cee(jn,ii) / m_ccc(jn,ii);
                m_cnw(jj,ii) = m_cnw(jj,ii) - m_cnn(jj,ii) * m_cww(jn,ii) / m_ccc(jn,ii);
                m_ccc(jj,ii) = m_ccc(jj,ii) - m_cnn(jj,ii) * m_css(jn,ii) / m_ccc(jn,ii);
                m_cnn(jj,ii) = -m_cnn(jj,ii) * m_cnn(jn,ii) / m_ccc(jn,ii);
            end
            if (js > 1)
                
                m_cse(jj,ii) = m_cse(jj,ii) - m_css(jj,ii) * m_cee(js,ii) / m_ccc(js,ii);
                m_csw(jj,ii) = m_csw(jj,ii) - m_css(jj,ii) * m_cww(js,ii) / m_ccc(js,ii);
                m_ccc(jj,ii) = m_ccc(jj,ii) - m_css(jj,ii) * m_cnn(js,ii) / m_ccc(js,ii);
                m_css(jj,ii) = -m_css(jj,ii) * m_css(js,ii) / m_ccc(js,ii);
            end
            if (ie <= m_Nx1 + 3)
                
                m_cne(jj,ii) = m_cne(jj,ii) - m_cee(jj,ii) * m_cnn(jj,ie) / m_ccc(jj,ie);
                m_cse(jj,ii) = m_cse(jj,ii) - m_cee(jj,ii) * m_css(jj,ie) / m_ccc(jj,ie);
                m_ccc(jj,ii) = m_ccc(jj,ii) - m_cee(jj,ii) * m_cww(jj,ie) / m_ccc(jj,ie);
                m_cee(jj,ii) = -m_cee(jj,ii) * m_cee(jj,ie) / m_ccc(jj,ie);
            end
            if (iw > 1)
                
                m_cnw(jj,ii) = m_cnw(jj,ii) - m_cww(jj,ii) * m_cnn(jj,iw) / m_ccc(jj,iw);
                m_csw(jj,ii) = m_csw(jj,ii) - m_cww(jj,ii) * m_css(jj,iw) / m_ccc(jj,iw);
                m_ccc(jj,ii) = m_ccc(jj,ii) - m_cww(jj,ii) * m_cee(jj,iw) / m_ccc(jj,iw);
                m_cww(jj,ii) = -m_cww(jj,ii) * m_cww(jj,iw) / m_ccc(jj,iw);
            end
        end
    end
    
    for jj = st+1:tw:m_Nx2+3
        for ii = st+1:tw:m_Nx1+3
            m_ccc(jj,ii) = m_ccc(jj,ii) +m_cnn(jj,ii) +m_cee(jj,ii) +m_css(jj,ii) +m_cww(jj,ii);
            
            %                 // Om geheugen te besparen gaan we cnw(jj,ii) opslaan in cnn(jj,ii);
            %                 //                                 cne(jj,ii) opslaan in cee(jj,ii);
            %                 //                                 csw(jj,ii) opslaan in cww(jj,ii);
            %         %                 //                                 cse(jj,ii) opslaan in css(jj,ii);
            %         m_cnn(jj,ii) = m_cnw(jj,ii);
            %         m_cee(jj,ii) = m_cne(jj,ii);
            %         m_cww(jj,ii) = m_csw(jj,ii);
            %         m_css(jj,ii) = m_cse(jj,ii);
            m_cnn(jj,ii) = 0;
            m_cee(jj,ii) = 0;
            m_cww(jj,ii) = 0;
            m_css(jj,ii) = 0;
            
        end
    end
    niv=niv+1;
    
    %
    %             // Uitgaande van het stencil
    %         //                          cnn
    %         //                     cnw      cne
    %         //                 cww      ccc      cee
    %         //                     csw      cse
    %         //                          css
    %         // elimineren we cnw en cne met de scheve 5-punts formule. Dit geeft
    %         //                 cnw      cnn      cne
    %         //                 cww      ccc      cee
    %         //                 csw      css      cse
    %         // We zijn nu weer terug de eerste 9-punts formule, maar op een
    %         // grover rooster.
    for jj = 1:tw:m_Nx2+3
        for ii = 1:tw:m_Nx1+3
            iw = ii - st;
            ie = ii + st;
            js = jj - st;
            jn = jj + st;
            if (jn <= m_Nx2 + 3 && ie <= m_Nx1 + 3)
                
                m_cnn(jj,ii) = m_cnn(jj,ii) - m_cne(jj,ii) * m_cnw(jn,ie) / m_ccc(jn,ie);
                m_cee(jj,ii) = m_cee(jj,ii) - m_cne(jj,ii) * m_cse(jn,ie) / m_ccc(jn,ie);
                m_ccc(jj,ii) = m_ccc(jj,ii) - m_cne(jj,ii) * m_csw(jn,ie) / m_ccc(jn,ie);
                m_cne(jj,ii) = -m_cne(jj,ii) * m_cne(jn,ie) / m_ccc(jn,ie);
            end
            if (iw > 1 && jn <= m_Nx2 + 3)
                
                m_cnn(jj,ii) = m_cnn(jj,ii) - m_cnw(jj,ii) * m_cne(jn,iw) / m_ccc(jn,iw);
                m_cww(jj,ii) = m_cww(jj,ii) - m_cnw(jj,ii) * m_csw(jn,iw) / m_ccc(jn,iw);
                m_ccc(jj,ii) = m_ccc(jj,ii) - m_cnw(jj,ii) * m_cse(jn,iw) / m_ccc(jn,iw);
                m_cnw(jj,ii) = -m_cnw(jj,ii) * m_cnw(jn,iw) / m_ccc(jn,iw);
            end
            if (js > 1 && iw > 1)
                
                m_css(jj,ii) = m_css(jj,ii) - m_csw(jj,ii) * m_cse(js,iw) / m_ccc(js,iw);
                m_cww(jj,ii) = m_cww(jj,ii) - m_csw(jj,ii) * m_cnw(js,iw) / m_ccc(js,iw);
                m_ccc(jj,ii) = m_ccc(jj,ii) - m_csw(jj,ii) * m_cne(js,iw) / m_ccc(js,iw);
                m_csw(jj,ii) = -m_csw(jj,ii) * m_csw(js,iw) / m_ccc(js,iw);
            end
            if (js > 1 && ie <= m_Nx1 + 3)
                
                m_css(jj,ii) = m_css(jj,ii) - m_cse(jj,ii) * m_csw(js,ie) / m_ccc(js,ie);
                m_cee(jj,ii) = m_cee(jj,ii) - m_cse(jj,ii) * m_cne(js,ie) / m_ccc(js,ie);
                m_ccc(jj,ii) = m_ccc(jj,ii) - m_cse(jj,ii) * m_cnw(js,ie) / m_ccc(js,ie);
                m_cse(jj,ii) = -m_cse(jj,ii) * m_cse(js,ie) / m_ccc(js,ie);
            end
        end
    end
    st = tw;
    tw = tw*2;
    
    
    
end

% for jj = 1:1:m_Nx2+3
%     for ii = 1:1:m_Nx1+3
%         
%         m_ccc(jj,ii) = 1.0 / m_ccc(jj,ii);
%         m_cnn(jj,ii) = m_cnn(jj,ii) * m_ccc(jj,ii);
%         m_css(jj,ii) = m_css(jj,ii) * m_ccc(jj,ii);
%         m_cee(jj,ii)= m_cee(jj,ii) * m_ccc(jj,ii);
%         m_cww(jj,ii) = m_cww(jj,ii) * m_ccc(jj,ii);
%         m_cnw(jj,ii) = m_cnw(jj,ii) * m_ccc(jj,ii);
%         m_cne(jj,ii) = m_cne(jj,ii) * m_ccc(jj,ii);
%         m_csw(jj,ii) = m_csw(jj,ii) * m_ccc(jj,ii);
%         m_cse(jj,ii) = m_cse(jj,ii) * m_ccc(jj,ii);
%         
%         
%     end
% end


end

